library(devtools)
library(obsval)
library(mvtnorm)
library(nnet)
library(MASS)
library(usethis)
library(haven)
library(ggplot2)
library(dplyr)
library(stargazer)
library(forecast)
library(tseries)


#load my data
myData <- read_dta("sop2021JLCdata.dta")



#####################
######## VFP ########
#####################

#vfp Negative Binomial
vfpNBreg <-  glm.nb(totalGrantedFJC ~ totalGrantedFJC_t1 + post1988 + diff_numCertPool + 
                      white + vfpDist + bindergrid1_t5 + diff_absAvgCircSCdist + diff_totalCertPet +
                      diff_pcthighestcoalition_t1 + diff_pctDisconnected_t1, data=myData)
vfpNBreg %>% summary()


#check residuals
# B-G: p-val less than .05 = independent
checkresiduals(vfpNBreg, 1, df = 55, plot = TRUE)
# Box: p-val greater than .05 = independent 
Box.test(resid(vfpNBreg),type="Ljung",lag=1)
adf.test(vfpNBreg$residuals, k=1)


#get pred values using obsval package
vfpPred <- obsval(totalGrantedFJC ~ totalGrantedFJC_t1 + post1988 + diff_numCertPool + 
                    white + vfpDist + bindergrid1_t5 + diff_absAvgCircSCdist + diff_totalCertPet +
                    diff_pcthighestcoalition_t1 + diff_pctDisconnected_t1, data=myData,
                  reg.model = "negbin",
                  n.draws = 1000,
                  effect.var = "vfpDist",
                  effect.vals = c(0, .07, .14, .21, .28, .35, .42),
                  #subsample.var = "post1988",
                  #subsample.val = 0,
                  verbose = TRUE)

#see results
vfpPred$model

#see effect of vfpconstraint
vfpPred$effect_sum

mean(vfpPred[["effect.preds"]])

vfpPred[["means"]]
vfpPred[["low.ci"]]
vfpPred[["high.ci"]]


vfpPredDat <- data.frame("range" = c(0, .07, .14, .21, .28, .35, .42),
                         "est" = c(133.8054, 130.9642, 128.2637, 125.5762, 122.9825, 120.6263, 118.0315),
                         "lb" = c(127.3521, 125.5589, 122.2831, 117.8809, 114.4403, 110.6767, 106.3195),
                         "ub" = c(140.1055, 137.0517, 134.3572, 132.7636, 131.3076, 130.7466, 129.3915))


vfpPredPlot<-ggplot(data=vfpPredDat, aes(x=vfpPredDat$range, y=vfpPredDat$est)) + 
  geom_line(size=3) +
  geom_errorbar(aes(ymin=vfpPredDat$lb, ymax=vfpPredDat$ub), width = .01, size=1) +
  #geom_line(aes(linetype=vfppremodel.df$div), size = 1.0) +
  ylab(label="Predicted Number of Cases on the Court's Docket\n") + 
  xlab(label="\nVeto-Filibuster Level of Constraint") +
  scale_linetype_manual(values = c(rep("solid"), rep("solid"))) +
  #ggtitle("Veto-Filibuster Pivot") +
  scale_x_continuous(limits = c(-.01, .43), breaks = c(0, .07, .14, .21, .28, .35, .42))+
  scale_y_continuous(limits = c(100, 150), breaks = c(100, 110, 120, 130, 140, 150))+
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5,face="bold", size=20)) +
  theme(axis.text=element_text(colour="black",size=20),
        axis.title=element_text(size=20)) +
  theme(panel.background = element_rect(colour = "black", size=1)) +
  theme(text = element_text(size=20)) +
  #theme(axis.title=element_text(face="bold")) +
  #theme(panel.grid.minor.x = element_blank()) +
  theme(panel.grid.major = element_line(colour = "gray67")) +
  theme(legend.title = element_blank()) 
vfpPredPlot



########################
######## MEDIAN ########
########################

#med Negative Binomial
medNBreg <-  glm.nb(totalGrantedFJC ~ totalGrantedFJC_t1 + post1988 + diff_numCertPool + 
                      white + medDist + bindergrid1_t5 + diff_absAvgCircSCdist + diff_totalCertPet +
                      diff_pcthighestcoalition_t1 + diff_pctDisconnected_t1, data=myData)
medNBreg %>% summary()

#check residuals
# B-G: p-val less than .05 = independent
checkresiduals(medNBreg, 1, df = 55, plot = TRUE)
# Box: p-val greater than .05 = independent 
Box.test(resid(medNBreg),type="Ljung",lag=1)
adf.test(medNBreg$residuals, k=1)


#get pred values using obsval package
medPred <- obsval(totalGrantedFJC ~ totalGrantedFJC_t1 + post1988 + diff_numCertPool + 
                    white + medDist + bindergrid1_t5 + diff_absAvgCircSCdist + diff_totalCertPet +
                    diff_pcthighestcoalition_t1 + diff_pctDisconnected_t1, data=myData,
                  reg.model = "negbin",
                  n.draws = 1000,
                  effect.var = "medDist",
                  effect.vals = c(0, .06, .12, .18, .24, .30, .36, .42, .48),
                  #subsample.var = "post1988",
                  #subsample.val = 0,
                  verbose = TRUE)


#see results
medPred$model

#see effect of medconstraint
medPred$effect_sum


medPred[["means"]]
medPred[["low.ci"]]
medPred[["high.ci"]]


medPredDat <- data.frame("range" = c(0, .06, .12, .18, .24, .30, .36, .42, .48),
                         "est" = c(128.9959, 127.1205, 125.5219, 123.7615, 122.2054, 120.5584, 118.9132, 117.3027, 115.7371),
                         "lb" = c(122.9981, 121.7335, 120.000, 117.9369, 115.7819, 113.1227, 110.7819, 107.5992,  105.735),
                         "ub" = c(135.3538, 133.0938, 131.2619, 129.4031, 128.7554, 127.9238, 127.1412, 126.3242, 125.9888))


medPredPlot<-ggplot(data=medPredDat, aes(x=medPredDat$range, y=medPredDat$est)) + 
  geom_line(size=1) +
  geom_rug(data=myData, aes(x=medDist, alpha = 1), size=1.25, sides="b", inherit.aes = FALSE, show.legend = FALSE) +
  geom_errorbar(aes(ymin=medPredDat$lb, ymax=medPredDat$ub), width = .01, size=1) +
  #geom_line(aes(linetype=vfppremodel.df$div), size = 1.0) +
  ylab(label="Predicted Number of Cases on the Court's Docket") + 
  xlab(label="Chamber Median Level of Constraint") +
  scale_linetype_manual(values = c(rep("solid"), rep("solid"))) +
  #ggtitle("Veto-Filibuster Pivot") +
  scale_x_continuous(limits = c(-.01, .49), breaks = c(0, .06, .12, .18, .24, .30, .36, .42, .48))+
  scale_y_continuous(limits = c(100, 140), breaks = c(100, 110, 120, 130, 140))+
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5,face="bold", size=20)) +
  theme(axis.text=element_text(colour="black",size=20),
        axis.title=element_text(size=20)) +
  theme(panel.background = element_rect(colour = "black", size=1)) +
  theme(text = element_text(size=20)) +
  #theme(axis.title=element_text(face="bold")) +
  #theme(panel.grid.minor.x = element_blank()) +
  theme(panel.grid.major = element_line(colour = "gray67")) +
  theme(legend.title = element_blank()) 
medPredPlot

ggsave(plot = medPredPlot, "/Users/elizabethlane/Dropbox/MPSA Revisions/2019 LaTex Files/1. JLC Revisions/LaTexFiles/medPredPlot.pdf")


###########################
######## JC MEDIAN ########
###########################

#####################
#judiciary committee median Negative Binomial
jcMedNBreg <-  glm.nb(totalGrantedFJC ~ totalGrantedFJC_t1 + post1988 + diff_numCertPool + 
                        white + jcMedDist + bindergrid1_t5 + diff_absAvgCircSCdist + diff_totalCertPet +
                        diff_pcthighestcoalition_t1 + diff_pctDisconnected_t1, data=myData)
jcMedNBreg %>% summary()


#check residuals
# B-G: p-val less than .05 = independent
checkresiduals(jcMedNBreg, 1, df = 55, plot = TRUE)
# Box: p-val greater than .05 = independent 
Box.test(resid(jcMedNBreg),type="Ljung",lag=1)
adf.test(jcMedNBreg$residuals, k=1)


#get pred values using obsval package
jcMedPred <- obsval(totalGrantedFJC ~ totalGrantedFJC_t1 + post1988 + diff_numCertPool + 
                      white + jcMedDist + bindergrid1_t5 + diff_absAvgCircSCdist + diff_totalCertPet +
                      diff_pcthighestcoalition_t1 + diff_pctDisconnected_t1, data=myData,
                    reg.model = "negbin",
                    n.draws = 1000,
                    effect.var = "jcMedDist",
                    effect.vals = c(0, .07, .14, .21, .28, .35, .42),
                    #subsample.var = "post1988",
                    #subsample.val = 0,
                    verbose = TRUE)
#see results
jcMedPred$model

#see effect of jcMedconstraint
jcMedPred$effect_sum


jcMedPred[["means"]]
jcMedPred[["low.ci"]]
jcMedPred[["high.ci"]]


jcMedPredDat <- data.frame("range" = c(0, .07, .14, .21, .28, .35, .42),
                           "est" = c(133.2588, 131.2104, 129.2695,  127.332, 125.4327, 123.6825, 121.7491),
                           "lb" = c(126.8475, 125.5237, 123.5411, 120.3877, 118.0097, 115.3216, 111.3386),
                           "ub" = c(139.5593, 137.4415, 135.3928, 134.4242, 132.7801, 132.3233, 132.1487))


jcMedPredPlot<-ggplot(data=jcMedPredDat, aes(x=jcMedPredDat$range, y=jcMedPredDat$est)) + 
  geom_line(size=3) +
  geom_errorbar(aes(ymin=jcMedPredDat$lb, ymax=jcMedPredDat$ub), width = .01, size=1) +
  #geom_line(aes(linetype=vfppremodel.df$div), size = 1.0) +
  ylab(label="Predicted Number of Cases on the Court's Docket\n") + 
  xlab(label="\nJudiciary Committee Median Level of Constraint") +
  scale_linetype_manual(values = c(rep("solid"), rep("solid"))) +
  #ggtitle("Veto-Filibuster Pivot") +
  scale_x_continuous(limits = c(-.01, .43), breaks = c(0, .07, .14, .21, .28, .35, .42))+
  scale_y_continuous(limits = c(100, 150), breaks = c(100, 110, 120, 130, 140, 150))+
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5,face="bold", size=20)) +
  theme(axis.text=element_text(colour="black",size=20),
        axis.title=element_text(size=20)) +
  theme(panel.background = element_rect(colour = "black", size=1)) +
  theme(text = element_text(size=20)) +
  #theme(axis.title=element_text(face="bold")) +
  #theme(panel.grid.minor.x = element_blank()) +
  theme(panel.grid.major = element_line(colour = "gray67")) +
  theme(legend.title = element_blank()) 
jcMedPredPlot





###################################
######## PARTY GATEKEEPING ########
###################################

#####################
#majority party Negative Binomial
majPartyNBreg <-  glm.nb(totalGrantedFJC ~ totalGrantedFJC_t1 + post1988 + diff_numCertPool + 
                           white + majPartyDist + bindergrid1_t5 + diff_absAvgCircSCdist + diff_totalCertPet +
                           diff_pcthighestcoalition_t1 + diff_pctDisconnected_t1, data=myData)
majPartyNBreg %>% summary()

#check residuals
# B-G: p-val less than .05 = independent
checkresiduals(majPartyNBreg, 1, df = 55, plot = TRUE)
# Box: p-val greater than .05 = independent 
Box.test(resid(majPartyNBreg),type="Ljung",lag=1)
adf.test(majPartyNBreg$residuals, k=1)


#get pred values using obsval package
majPartyPred <- obsval(totalGrantedFJC ~ totalGrantedFJC_t1 + post1988 + diff_numCertPool + 
                         white + majPartyDist + bindergrid1_t5 + diff_absAvgCircSCdist + diff_totalCertPet +
                         diff_pcthighestcoalition_t1 + diff_pctDisconnected_t1, data=myData,
                       reg.model = "negbin",
                       n.draws = 1000,
                       effect.var = "majPartyDist",
                       effect.vals = c(0, .07, .14, .21, .28, .35, .42),
                       #subsample.var = "post1988",
                       #subsample.val = 0,
                       verbose = TRUE)

#see results
majPartyPred$model

#see effect of majPartyconstraint
majPartyPred$effect_sum


jcMedPred[["means"]]
jcMedPred[["low.ci"]]
jcMedPred[["high.ci"]]

majPartyPredDat <- data.frame("range" = c(0, .07, .14, .21, .28, .35, .42),
                              "est" = c(133.2588, 131.2104, 129.2695,  127.332, 125.4327, 123.6825, 121.7491),
                              "lb" = c(126.8475, 125.5237, 123.5411, 120.3877, 118.0097, 115.3216, 111.3386),
                              "ub" = c(139.5593, 137.4415, 135.3928, 134.4242, 132.7801, 132.3233, 132.1487))


majPartyPredPlot<-ggplot(data=majPartyPredDat, aes(x=majPartyPredDat$range, y=majPartyPredDat$est)) + 
  geom_line(size=3) +
  geom_errorbar(aes(ymin=majPartyPredDat$lb, ymax=majPartyPredDat$ub), width = .01, size=1) +
  #geom_line(aes(linetype=vfppremodel.df$div), size = 1.0) +
  ylab(label="Predicted Number of Cases on the Court's Docket\n") + 
  xlab(label="\nMarjority Party Median Level of Constraint") +
  scale_linetype_manual(values = c(rep("solid"), rep("solid"))) +
  #ggtitle("Veto-Filibuster Pivot") +
  scale_x_continuous(limits = c(-.01, .43), breaks = c(0, .07, .14, .21, .28, .35, .42))+
  scale_y_continuous(limits = c(100, 150), breaks = c(100, 110, 120, 130, 140, 150))+
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5,face="bold", size=20)) +
  theme(axis.text=element_text(colour="black",size=20),
        axis.title=element_text(size=20)) +
  theme(panel.background = element_rect(colour = "black", size=1)) +
  theme(text = element_text(size=20)) +
  #theme(axis.title=element_text(face="bold")) +
  #theme(panel.grid.minor.x = element_blank()) +
  theme(panel.grid.major = element_line(colour = "gray67")) +
  theme(legend.title = element_blank()) 
majPartyPredPlot




##################################
######## PRINT OUT TABLES ########
##################################

stargazer(vfpNBreg, medNBreg, jcMedNBreg, majPartyNBreg, title="Negative Binomial Regression Models of Supreme Court Caseload",
          align=TRUE, dep.var.labels=c("Veto-Fillibuster","Chamber Median", "Committee Gatekeeping", "Party Gatekeeping"),
          covariate.labels=c("Total Granted t-1", "Post 1988", "Diff Cert. Petitions",
                             "Diff Ideological Range", "Justice White",  "vfp", "med", "jcmed", "majparty", "bindergrid1_t5", 
                             "diff_absAvgCircSCdist",  "diff_totalCertPet"), omit.stat=c("LL","ser","f"), no.space=TRUE)





##################################
######## APPENDIX TABLE 2 ########
##################################
adf.test(myData$vfpDist)
kpss.test(myData$vfpDist)

adf.test(myData$medDist)
kpss.test(myData$medDist)

adf.test(myData$jcMedDist)
kpss.test(myData$jcMedDist)

adf.test(myData$majPartyDist)
kpss.test(myData$majPartyDist)

adf.test(myData$pcthighestcoalition_t1)
kpss.test(myData$pcthighestcoalition_t1)
adf.test(myData$diff_pcthighestcoalition_t1)
kpss.test(myData$diff_pcthighestcoalition_t1)

adf.test(myData$pctDisconnected_t1)
kpss.test(myData$pctDisconnected_t1)
adf.test(myData$diff_pctDisconnected_t1)
kpss.test(myData$diff_pctDisconnected_t1)

adf.test(myData$numCertPool)
kpss.test(myData$numCertPool)
adf.test(myData$diff_numCertPool)
kpss.test(myData$diff_numCertPool)

adf.test(myData$absAvgCircSCdist)
kpss.test(myData$absAvgCircSCdist)
adf.test(myData$diff_absAvgCircSCdist)
kpss.test(myData$diff_absAvgCircSCdist)

adf.test(myData$totalCertPet)
kpss.test(myData$totalCertPet)
adf.test(myData$diff_totalCertPet)
kpss.test(myData$diff_totalCertPet)




